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Abstract: Two methods for the determination of extrema timings and 
their uncertainties appropriate for the analysis of time series of variable 
stars using matrix calculus are presented. The method I is suitable 
for determination of times of extrema of non-periodical variables or 
objects, whose light curves vary. The method II is apt for O-C analyses 
of objects whose light curves are more or less repeating. 
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1. Linear regression using matrix calculus 

The determination of times of extrema and its uncertainty of light vari- 
ations belongs to the most common problems in the astrophysics of variable 
stars. The following instructive text develops the basic ideas outlined in jl|) 
Mikulasek et al. (2006). 

1.1. The data 

Let we have n observations {y^} creating the vector Y, Y = {yi} done in 
moments {ij} creating the vector t, t = {ti}, both of the orders of n x 1. It 
is demanding to adapt the time scale so that t =>• t — t. 

The weight Wi expresses the quality of an i-th measurement; w = {wi}. 
If we know the inner uncertainty (error) of the i-th observation 5yi, we can 
determine the weight of it as follows: Wi ~ (5yi)~ 2 . If we assume that the 
quality of all measurements are more or less equal or if we know about the 
quality of measurements little or nothing, it would be honest to put Wi = 1, 
w = E(n, 1). 

The most of relations look simpler if we normalise weights so it is valid 
W = 1, then Ya=1 w i = n - Frequently, we use instead of the w square matrix 
(n x n) W, W = diag(w). 



1.2. Linear regression 

The observed relationship between the dependent variable (inaccurately me- 
asured quantity - mostly magnitude, radial velocity, temperature) y and the 
independent variable (precisely measured quantity - typically time) t can be 
fit by an appropriate model function F(t). The model function is determi- 
ned by g free parameters (3j that create a column vector j3 = [Pi, 02, ...(3 g ] T . 
The upper index denotes transposing the matrix. If the model function 
F(t) can be expressed as a linear combination of g different functions oi 
time fk(t) (we speak here about the linear model function). Then 



f(t) = ih(t), / 2 (t), ■..,/*(*)], F(t, /3) = ]T Pkfk(t) = f(i)/3- (1) 

fc=i 

Let introduce the matrix X of rank nx g and the column vector Y p , (nxl): 



X = 



/ f(*i) \ 
V f(«n) y 



Y P = 



^(* 2 ) 



= X/3. 



(2) 



As the objective measure of the success rate of the fit for an a (3 is used 
usually the sum of weighted squares of deflection of observed values and 
predicted ones S((3): 



S(/3) = (Y- Y P ) T W(Y- Y p ) = Y T WY-2/3 T U + /3 T V/3, 
where 

U = X T WY; V = X T WX; H = V \ 



(3) 



(4) 



The least square method (LSM) considers the fit by the model function 
F(t, f3) as the best one if the sum R = S(/3 = b) is minimal. In the case oi 
the linear model function F(t,(3) we obtain for b and R: 



dS_ 



= = -2U + 2Vb, 



b = HU; 



R = Y T WY-b T U. 



(5) 



1.3. The example 



Standardly used linear regression models are polynomials or trigonometric 
functions of arbitrary orders. As an example we select the parabolic model 



- the simplest model we can use for fit of the real function in the form: 
F(t) = ft i 2 + ft t + ft, f (*) = [t 2 , t,l],X= M } {!}]. 



L^. Standard deviation. Uncertainties 
The standard deviation cr can be estimated using relation 




n-g' 



(6) 



The components of the column vector Sh used to be considered as a rigou- 
rous estimate of the uncertainty of the particular parameters. Unfortuna- 
tely, they have this meaning only exceptionally, nevertheless it is sometimes 
required by referees. Contrary, very valuable is the following estimate of the 
uncertainty of the model predictions SF(t) 



Method I of the determination of the extrema times consists of several steps. 

• Plot the time series and select the appropriate linear model function 
which can fit the observed light curve sufficiently well with the mini- 
mum of free parameters. 

• Fit the observed light curve by the model function and check whether 
the dependence of deflections of observed and predicted light curves 
show only random scatter or some trends. In the latter case improve 
your model function 

• Find the extreme time on the fitted light curve and calculate its un- 



2.1. When the minimum /maximum occurs? 

The moments of the fitted model function extrema occurs if it is fulfilled 
the condition that time derivative in that t e equals to zero, especially 



5b = a v /diag(H); SF(t) = a ^f(i)Hf(t) T . 



(7) 



2 Method I 



certainty. 



= y p (t c ) = F(t e ) = 



d(f(t c )b) 



= f(*e) 



b; ¥(t) = [f 1 (t)J 2 (t),...f g (t)]. (8) 



dt 



There are many techniques how to find roots of the derivative of the fitted 
function F(t). The uncertainty of the determination of the time derivative 
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Obrazek 1: The simulated light curve of a variable star with a minimum. Inner 
accuracy of individual measurements are denoted by gray error bars. The fitted 
parabola is marked by black dots with error bars corresponding to the uncertainty 
of predicted values. 



at the time t, 5y p (t) can be estimated by the relation: 



5y p (t) = ayJf(t)Hf(t) T - 



0) 



For the estimate of the uncertainty of the time of extreme St e we need yet 
the second time derivative in of the fitted curve in the time of extreme 

y P (Q- 

y p (Q = f(t e )b; m = [f 1 (t)J 2 (t),...f g (t)] ] 5t e = ^^. (10) 



3 Method II 



This method is again an application of the weighted least square method 
this once with the model function defined so that time or times of extrema 
are free parameters which are found and iteratively improved. Uncertainties 



of these parameters are then simply uncertainties of determined extrema 
moments. 

The method is applied namely for periodic or nearly periodic variable 
stars as pulsating stars, eclipsing binaries or rotating stars. The disadvan- 
tage of the method consists in the fact that from principle reasons we cannot 
applied here linear regression a we have to used approaches established for 
solution of non-linear LSM. 

3.1. Example - transformation of the parabolic model 

In section 1.3 we introduced the linear model expression with three para- 
meters /?, F(t) = j3\ t 2 + /?2 t + (3%. This parabolic model function can be 
rewritten in the non-linear form: 

F(t, a) = a 2 {t - ai) 2 + a 3 , (11) 

where the parameter a\ = t m i n corresponds to the time of the minimum oi 
quadratic model function, a\ = —fa/^Pi, «2 = Pi, oc^ = Ps — /3f/4/3i. It 
is apparent that the model F(t, a) as it is expressed in Eq.QT] is no more 
linear, nevertheless it is in the form demanded by the method II. 

3.2. Linearisation of non-linear model functions 

The solution of non linear regression is not straightforward and immediate 
than in the case of linear regression. Nevertheless, having a satisfactorily 
an initial estimate of the solution clq we can linearise the non-linear mo- 
del function and delight in all advantages yielding by linear models. The 
linearisation can be done by the Taylor decomposition of the first order 




The function is now linear in respect to new set parameters Aa, functions 
of the linearised model are f(i) = [^^, ■■■^ L ]- 

Now we can create matrices X (see Eq.[2]), Y, W and calculate a solution 
for the vector Aa. 

V = X T WX; U = X T WX; H = V -1 ; Aa = HU. (13) 

This correction we add to the initial estimate ao and we obtained the further, 
improved estimate for the vector of parameters ai = ao + Aa. The new value 
of parameter set we can repeat the whole procedure. The iterative process 
diminish the value of the correcting vector Aa as rule very effectively and 
after several steps we obtain the final result. 



You need not to iterate all parameters, some of them (the linear ones) 
can be calculate directly - see the following example. 

3.3. Example - linearisation of the parabolic model 

Now we can linearise our quadratic model according to the Eq.[12l assuming 
the initial estimate of parameters ao 

F= a 02 (t - a i) 2 + a 3 + Aoi 2 a 02 (t - a i) + Aa 2 (t - a i) 2 + Aa 3 . (14) 

The function is in parameters a 2 , a 3 linear. It means we can calculate the 
first two parameters directly, and a\ iteratively: 

F(t, a) = Aa\ 2a 2 (t — a i) + a 2 (t - a i) 2 + a 3 . (15) 



( 2a 2 (h - a i) (h 
x _ 2a 2 (t2-aoi) (h 

\ 2a 2 (t n - a i) (t n 

Knowing parameters a 2 , 03 (repeating light variations) we should fix them. 
3-4- Estimation of the uncertainty of extrema timings 

The uncertainties of free parameters including the extremum timing are 
given by 

AY = Y — F(t, a); R = AY T W AY; 5a = o-^/diag(H). (17) 

Anyhow the resulting sets of correcting parameters are practically 'pure 
zero', their uncertainties differ from zero a correspond to uncertainties oi 
the particular parameters. It enables to determine the reliable estimate for 
errors in extrema time determinations. 



am) 2 1 \ 
aoi? 1 

am) 2 1 J 



(16) 



4 Discussion and conclusion 

It can be also proven that in the case of identical model the both methods 
are interchangeable and yield the identical results. 
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